function output = estimate(pars,M_full,W_full,Nu,verbose)     
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Vector of fixed parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rho                                     = NaN(numel(Nu.fixed_pars),1);
for p=1:numel(Nu.fixed_pars)
        Rho(p)                                  = pars.(Nu.fixed_pars{p});
end 

%%%%%%%%%%%%%%
% Constraints
%%%%%%%%%%%%%%
A                                       = Nu.A;
b                                       = Nu.b;
lb                                      = Nu.lb;
ub                                      = Nu.ub; 

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Grid of starting points
%%%%%%%%%%%%%%%%%%%%%%%%%%
if isequal(Nu.targeted_pars,{'S','M_e','C'})
S_vec                                   = linspace(lb(1),ub(1),Nu.Nvec)';
C_vec                                   = linspace(lb(3)+Nu.rh*(ub(3)-lb(3)),ub(3)-Nu.rh*(ub(3)-lb(3)),Nu.Nvec)';
[S_grid,C_grid,M_e_grid]                = deal(NaN(Nu.Nvec^3,1));
for n=1:Nu.Nvec
        S_grid( ((n-1)*(Nu.Nvec^2)+1):(n*(Nu.Nvec^2)) ) = S_vec(n);
        for m=1:Nu.Nvec
                C_grid( ((n-1)*(Nu.Nvec^2))+(((m-1)*(Nu.Nvec)+1):(m*(Nu.Nvec))) )       = C_vec(m);
                ub_M                                                                    = min(ub(2),pars.M_c + pars.xi*sqrt(Nu.N/pars.theta) - C_vec(m));
                if ub_M < lb(2)
                        save('runs/evaluate_call.mat');
                        error('Error: problem with lb and ub in grid \n');
                end
                M_e_grid( ((n-1)*(Nu.Nvec^2))+(((m-1)*(Nu.Nvec)+1):(m*(Nu.Nvec))) )     = linspace(lb(2),ub_M,Nu.Nvec)';
        end
end
Dummy                                   = [S_grid,M_e_grid,C_grid];
if any(A*(Dummy') > b)||any( any( Dummy > repmat(ub',[Nu.Nvec^3 1]) ) )||any( any( Dummy < repmat(lb',[Nu.Nvec^3 1]) ) )
        save('runs/estimate_call.mat');
        error('Error: point violates the local constraints \n');
end
Dummy                                   = uniquetol(Dummy,'ByRows',true);
Dummy                                   = sortrows(Dummy,[1 3 2],'ascend');
S_grid                                  = Dummy(:,1);
M_e_grid                                = Dummy(:,2);
C_grid                                  = Dummy(:,3);
Nevals                                  = numel(S_grid);
else
    error('No grid of starting points for this set of parameters');
end

%%%%%%%%%%%%%%%%%%%%%%
% Objective function
%%%%%%%%%%%%%%%%%%%%%%
fun                                     = @(x) SMM_obj( x, Rho, Nu, M_full, W_full, Nu.log_targeted_moments, 0 );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Evaluate at potential starting points
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ObjStart                                = NaN(Nevals,1);
XiStart                                 = NaN(Nevals,3);
parfor n=1:Nevals
        dummyXi                                 = [S_grid(n),M_e_grid(n),C_grid(n)];
        dummyObj                                = fun(dummyXi);
        ObjStart(n)                             = dummyObj;
        XiStart(n,:)                            = dummyXi;
end
AllStart                                = [(1:Nevals)',ObjStart,XiStart];
AllStart                                = sortrows(AllStart,2);
StartingPoints                          = AllStart(1:Nu.Nstart,2:end);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Optimize at best starting points
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ObjStart                                = StartingPoints(:,1);
XiStart                                 = StartingPoints(:,2:end)';
[ObjHat,Fevals,Flag]                    = deal(NaN(size(ObjStart)));
XiHat                                   = NaN(size(XiStart));
parfor n=1:Nu.Nstart
        Xi_0                            = XiStart(:,n);
        [Xi,Obj,loc_flag,other]         = patternsearch(fun,Xi_0,A,b,[],[],lb,ub,[],Nu.options_patternsearch);
        XiHat(:,n)                      = Xi;
        ObjHat(n)                       = Obj;
        Fevals(n)                       = other.funccount;
        Flag(n)                         = loc_flag;
end
AllHat                                  = [(1:Nu.Nstart)',ObjHat,XiHat',Fevals,Flag];
AllHat                                  = sortrows(AllHat,2);
PointFinal                              = AllHat(1,:);
ObjFinal                                = PointFinal(2);
XiFinal                                 = PointFinal(3:5);
FevalsFinal                             = PointFinal(6);
FlagFinal                               = PointFinal(7);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute objective and moments at final point
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dummy                                   = produce_moments(XiFinal,Rho,Nu);
OutputFinal.M                           = Dummy.M;
OutputFinal.X_monthly                   = Dummy.X_monthly;
OutputFinal.M_monthly                   = Dummy.M_monthly;
ObjFinalAllMoms                         = SMM_obj( XiFinal, Rho, Nu, M_full, W_full, true(size(Nu.log_targeted_moments)), 0 );
XiFinalAllPars                          = [pars.xi,pars.sigma,pars.k,XiFinal(1),XiFinal(2),XiFinal(3)]';

%%%%%%%%%%%%%%%
% Print output
%%%%%%%%%%%%%%%
if verbose 
    Table6                                  = table(num2str(XiFinalAllPars,'%.5f'),...
                                                    strvcat(num2str(ObjFinal,'%.4f'),repmat(' ',[numel(XiFinalAllPars)-1 1])),...
                                                    strvcat(num2str(ObjFinalAllMoms,'%.4f'),repmat(' ',[numel(XiFinalAllPars)-1 1])),...
                                                    strvcat(num2str(FlagFinal,'%.0f'),repmat(' ',[numel(XiFinalAllPars)-1 1])),...
                                                    strvcat(num2str(FevalsFinal,'%.0f'),repmat(' ',[numel(XiFinalAllPars)-1 1])));
    Table6.Properties.VariableNames         = ["Point estimate" "Obj. (targ.)" "Obj. (all)" "Flag" "Fevals"];
    Table6.Properties.RowNames              = {'xi','sigma','k','S','M_e','C'};
    disp(Table6)
    Table7                                  = table(num2str(M_full,'%.3f'),...
                                                    num2str(mean(OutputFinal.M,2),'%.3f'));
    Table7.Properties.VariableNames         = ["Emp. value" "Sim. value"];
    Table7.Properties.RowNames              = Nu.all_moments;
    disp(Table7)
end

%%%%%%%%%
% Output
%%%%%%%%%

output.XiFinalAllPars                       = XiFinalAllPars;
output.ObjFinal                             = ObjFinal;
output.ObjFinalAllMoms                      = ObjFinalAllMoms;
output.FlagFinal                            = FlagFinal;
output.FevalsFinal                          = FevalsFinal;
output.OutputFinal                          = OutputFinal;

end